TR-3tia 


l£  C 


TECHNICAL  REPORT  BRL-TR-2915 


BRL 


1938  -  Serving  the  Army  for  Fifty  Years  —  1988 


NUMERICAL  SIMULATION  OF  UNSTEADY 
INCOMPRESSIBLE  FLOW  IN  A 
PARTIALLY-FILLED  ROTATING  CYLINDER 


MICHAEL  J.  NUSCA 


DTIC 

e:lecte| 

JUN  1  7  19881 


“b 


JUNE  1988 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


U.S.  ARMY  LABORATORY  COMMAND 


BALLISTIC  RESEARCH  LABORATORY 


ABERDEEN  PROVING  GROUND,  MARYLAND 


38  6  IS  0  5  0' 


DESTRUCTION  NOTICE 


Destroy  this  report  when  it  is  no  longer  needed.  DO  NOT  return  it  to  the 
originator . 


Additional  copies  of  this  report  may  be  obtained  from  the  National  Technical 
Infonretion  Service,  U.S.  Department  of  Cattnerce,  Springfield,  VA  22161. 


The  findings  of  this  report  are  not  to  be  construed  as  an  official  Department 
"of  the  Army  position,  unless  so  designated  by  other  'authorized  documents. 


The  use  of  trade  names  or  manufacturers'  names  in  this  report  does  not  con¬ 
stitute  indorsement  of  any  carmercial "  product .  .  ■■■-. 


why.  y 


8c.  ADDRESS  (City,  State,  and  ZIP  Code) 

Aberdeen  Proving  Ground,  MD  21005-5066 


1 1  TITLE  (Include  Security  Classification) 

Numerical  Simulation  of  Unsteady  Incompressible  Flow  in  a  Partially-Filled  Rotating  Cylir 


|  10  SOURCE  OF  FUNDING  NUMBERS 

PROGRAM 
ELEMENT  NO. 

62618A 

PROJECT 

N0  1L1 
62618AH80 

TASK 

NO 

WORK  UNIT 
ACCESSION  N 

GROUP 


SUB-GROUP 


Finite  Difference  •  Liquid  Moment  ' 

Incompressible  Flow  '  Unsteady  Flow 

LiquidJFilled  Proiectile  .  Rotating  Liquids,  h  '  . 


19  ABSTRAiQT  ( Continue  on  reverse  if  necessary  and  identify  by  block  number) 

The  liquid  flowfield  in  a  full  or  partially-filled  right  circular  cylinder  in  rapid 
axial  rotation  is  investigated  numerically.  The  governing  equations  are  the  axisymmetric i 
unsteady,  viscous,  incompressible  Navier-Stokes  equations.  These  equations  are  written  it 
stream  function-vorticity  form  for  a  cylindrical  coordinate  system  in  a  nonrotating  refer¬ 
ence  frame.  The  governing  equations  are  discretized  using  second-order  finite-difference* 
for  time  and  space  on  a  nonuniform  grid  employing  logarithmic  stretching  in  regions  where 
high  flow  gradients  are  anticipated.  Time  dependent  solutions  for  Reynolds  numbers  betwei 
1,000  and  100,000  have  been  obtained  using  a  Gauss-Seidel  relaxation  procedure.  For  partJ 
ly  filled  cases  the  free  surface  is  assumed  to  be  cylindrical  and  located  at  a  constant 
radius  from  the  axis  of  spin.  Numerical  solutions  for  full  cylinders  are  consistent  with 
previous  solutions  and  experimental  data.  Numerical  solutions  for  a  partially-filled  cyli 
der  are  consistent  with  experimental  data  for  a  liquid  centrifuge  except  at  the  free  surf) 

Computations  of  the  roll  moment  exerted  on  the  cylinder  by  the  contained  liquid  shows  a 
smaller  moment  for  the  partially-filled  compared  with  the  full  cylinder  results.  !t  UU-r> 


20  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT 
□  UNCLASSIFIED/UNLIMITED  Q  SAME  AS  RPT  □  DTIC  USERS 


22«  NAME  OF  RESPONSIBLE  INDIVIDUAL 
Ml  rhap)  .1 .  Ntisca 


21  ABSTRACT  SECURITY  CLASSIFICATION 
UNCLASSIFIED 


22b  TELEPHONE  (Include  Area  Code)  22c.  OFFICE  SYMBOL 
(30D-278-2057  SLCBR-LF-a 


Acknowledgement 

The  following  individuals  have  contributed  to  the  completion  of  this  work: 
Dr.  William  P.  D’Amico,  Dr.  Gene  Cooper,  and  Mr.  Nathan  Gerber  of  the 
Launch  and  Flight  Division,  U.S.  Army  Ballistic  Research  Laboratory,  Aberdeen 
Proving  Ground,  Maryland.  Dr.  Sukumar  Chakravarthy  of  Rockwell  Interna¬ 
tional  Science  Center,  Thousand  Oaks,  California,  provided  technical  guidance  in 
the  implementation  of  the  numerical  solution  scheme  for  the  equations  of  motion. 
This  work  was  completed  in  partial  fulfillment  of  a  Master  of  Science  degree  at 
the  University  of  Maryland;  Dr.  John  D.  Anderson,  advisor. 


Accesion  For 

_ _ _ _ _ 

NT'S  CRA&I 
OTiC  TAB 

U'  Jl.tlQ  i  CC'.'J 

!  J.iMi'iCc'ti-.ji; 


I  <:U"/ 


!  4  ’  Av--  :  (. 

1  ■  ■ 


\A'i 


i  l 


Cory 

JIWltTK 


Table  of  Contents 


Page 


ACKNOWLEDGMENTS .  iii 

LIST  OF  FIGURES .  vii 

I.  INTRODUCTION .  1 

II.  FREE  SURFACE  CONSTRAINTS .  2 

III.  GOVERNING  EQUATIONS .  4 

IV.  BOUNDARY  CONDITIONS .  5 

V.  NUMERICAL  SOLUTION  SCHEME .  8 

VI.  COMPARISONS  WITH  PREVIOUS  WORK .  10 

VII.  RESULTS  AND  DISCUSSION .  13 

VIII.  CONCLUSIONS .  15 

REFERENCES .  43 

LIST  OF  SYMBOLS .  45 

DISTRIBUTION  LIST  .  47 


l\ 

<■" 

> 

A 


List  of  Figures 

Figure  Page 


1  Difference  in  free  surface  radial  location  at  the  top  wall  and  bottom  wall  of 

the  cylinder  as  a  function  of  Froude  number,  95%  fill .  16 

2  Cylinder  geometry  and  computational  domain .  17 

3  Finite-difference  grid  (31  x  31)  with  logarithmic  stretching  (d  =  b  =  0.05) 

for  full  cylinder  ( c/a  =  1.0) .  18 

4  Finite-difference  grid  (31  x  31)  with  logarithmic  stretching  (d  =  b  =  0.05) 

for  70%  filled  cylinder  ( c/a  =  1.0)  .  18 

5  Azimuthal  velocity  profiles  at  cylinder  midplane  for  a  =  1,  Re  =  9741.6  .  .  19 

6  Azimuthal  velocity  profiles  at  cylinder  midplane  for  a  =  1.515,  Re  =  3076  .  20 

7  Azimuthal  velocity  profile  near  cylinder  endwall  ( z  =  0.092)  for  a  =  1,  Re 

=  9741.6  .  .  .  .* .  21 

8  Radial  velocity  profile  in  Ekman  layer  at  r  =  0.76  for  a  =  1,  Re  =  9741.6 

and  t  =  4.5  sec .  22 

9  Inertial  oscillations  during  spin-up  from  a  previous  state  of  rigid  body  rota¬ 
tion  at  r  =  0.25,  z  =  0.3182 .  23 

10  Inertial  oscillations  during  spin-up  from  a  previous  state  of  rigid  body  rota¬ 
tion  at  r  =  0.25,  z  =  0.0790  .  24 

11  Inertial  oscillations  during  spin-up  from  a  previous  state  of  rigid  body  rota¬ 
tion  at  r  =  0.75,  z  =  0.3182 .  25 


£ 


r  a 
« 

p 

l 

8 


" » 

8 

i 

J 


12  Inertial  oscillations  during  spin-up  from  a  previous  state  of  rigid  body  rota¬ 
tion  at  r  =  0.75,  z  =  0.0790  .  26 

13  LDV  measurement  stations  for  liquid  centrifuge  experiment  (Ref.  9);  1:  z 

=  0.497,  2:  z  =  0.9,  3:  z  =  1.304,  4:  z  =  1.706  . . .  27 

14  Finite-difference  grid  (30  x  61)  for  Re  =  93518,  a  =  1.028,  36.7%  filled. 

(A z)endwatu  =  0.0018,  (A r)sidewa,i  =  (A r)/ree  surface  =  0.0014 .  27 

15  Axial  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  z  =  0.497  .  .  28 

16  Azimuthal  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  z  =  0.497  29 

17  Axial  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  z  =  0.9  ...  30 

18  Azimuthal  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  z  =  0.9  31 

19  Axial  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  z  =  1.304  .  .  32 

20  Azimuthal  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  z  =  1.304  33 


List  of  Figures  (Continued) 

Figure  Page 

21  Axial  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  z  =  1.706  .  .  34 

22  Azimuthal  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  2  =  1.706  35 


23  Finite-difference  grid  (21  x  41)  for  Re  =  1000,  a  —  3,  full  cylinder .  36 

24  Finite-difference  grid  (21  x  41)  for  Re  =  1000,  a  =  3,  95%  filled  cylinder  .  .  36 

25  Toted  liquid  roll  moment  coefficient  vs  time  for  Re  =  1000,  a  =  3,  100%  and 

95%  filled  cylinders .  37 

26  Sidewall  component  of  liquid  roll  moment  coefficient  vs  time  for  Re  =  1000, 

a  =  3,  100%  and  95%  filled  cylinders .  38 

27  Endwall  component  of  liquid  roll  moment  coefficient  vs  time  for  Re  =  1000, 

a  =  3,  100%  and  95%  filled  cylinders .  39 

28  Azimuthal  velocity  profiles  for  2  =  0.009  and  3.0  for  Re  =  1000,  c*  =  3,  t  = 

0.00025  sec,  100%  and  95%  filled  cylinders .  40 

29  Relative  difference  in  total  liquid  roll  moment  coefficient  for  /  =  1  and  /  = 

0.95  vs  time  for  Re  —  1000,  a  =  3 .  41 


I.  INTRODUCTION 


The  objective  of  this  work  was  to  investigate  the  axisymmetric,  unsteady,  viscous, 
incompressible  liquid  flow  in  a  rapidly  rotating  right  circular  cylinder.  Understanding  of 
this  internal  flow  is  necessary  in  the  design  of  gun-launched  projectiles  which  carry  liquid 
payloads.  Liquid-filled  projectiles  exhibit  flight  instabilities  that  are  not  observed  for  the 
same  vehicle  with  a  rigid  payload.  These  instabilities  have  been  described  by  Stewartson.1 
The  Stewartson  principle  states  that  resonance  between  the  inertial  wave  frequency  in  the 
rotating  liquid  and  the  nutation  frequency  of  the  projectile  can  cause  projectile  instability. 
A  numerical  procedure  for  solving  the  unsteady  Navier-Stokes  equations  was  developed 
for  the  present  study.  From  a  computational  viewpoint,  this  problem  is  instructive  since 
it  represents  a  class  of  internal  flow  problems  for  which  computational  experiments  can 
reveal  details  of  the  flow  that  cannot  be  easily  measured  experimentally.  This  study  of 
axisymmetric  free  surface  flow  represents  the  initial  stage  of  an  investigation  of  the  three- 
dimensional  problem. 

Theoretical  and  experimental  work  examining  the  flow  within  a  rotating  cylinder  has 
concentrated  on  full  cylinders.2  Experimental  evidence3  indicates  that  the  effect  of  a  free 
surface  on  enclosed  rotating  liquids  is  important.  Particular  interests  include  shear  layers 
that  may  form  on  the  liquid-air  interface,  the  role  of  liquid  viscosity  and  interior  flow 
circulation  in  restoring  the  disturbed  fluid  to  a  state  of  rigid  body  rotation,  the  pressure 
and  shear  forces  exerted  on  the  cylinder  walls  by  the  liquid,  and  the  resulting  liquid  roll 
moment  on  the  container. 

Computational  work  for  this  problem  has  taken  several  forms.  Kitchens4  adapted  a 
predictor-corrector,  multiple-iteration  scheme,  which  -was  developed  by  Rubin  and  Lin5  for 
steady  flows,  combined  with  a  Gauss-Seidel  iteration  scheme  to  solve  the  time-dependent 
flow  in  a  cylinder  spinning  up  from  rest.  Chakravarthy6  composed  two  time  accurate 
schemes  to  compute  the  time-dependent  flow  for  the  spin-up  problem.  In  one  scheme  the 
finite-differenced  equations  in  primitive  variable  form  were  solved  using  an  implicit  approxi¬ 
mate  factorization  procedure.  The  second  scheme  solved  the  equations  in  stream  function- 
vorticity  form  using  a  Gauss-Seidel  relaxation  procedure.  Kitchens  and  Chakravarthy 
produced  comparable  results  for  filled  cylinders.  Homicz  and  Gerber7  extended  a  model 
by  Goller  and  Ranov8  and  produced  a  finite-difference  scheme  for  spin-up  from  rest  in  a 
partially-filled  rotating  cylinder.  This  analysis  assumed  columnar  flow  (i.e.  flow  velocities 
and  pressure  are  independent  of  the  axial  coordinate)  and  thus  the  azimuthal  momentum 
equation  was  the  sole  equation  considered.  This  equation  was  solved  using  an  implicit 
Crank-Nicolson  scheme. 

Shadday9  conducted  an  experimental  and  numerical  investigation  of  the  internal  flow- 
field  of  a  liquid  centrifuge.  He  considered  a  partially-filled  right  circular  cylinder  in  rapid 
rotation  with  the  internal  flowfield  induced  by  the  differential  rotation  of  one  endwall 
with  respect  to  the  rest  of  the  container.  Shadday  gathered  axial  and  azimuthal  velocity 
data  along  radii  of  the  cylinder  using  a  laser-Doppler  velocimeter  (LDV).  In  addition,  a 
numerical  scheme  was  formulated  in  primitive  variables  and  used  a  modified  marker-and- 
cell  (MAC)  implicit  scheme.  This  scheme  included  an  explicit  treatment  of  the  governing 
equations  with  the  exception  of  the  Coriolis  acceleration  terms.  The  numerical  study  used 


a  coarse  grid  in  the  radial  direction.  An  inadequate  correlation  between  numerical  and 
experimental  data  at  the  liquid  free  surface  was  noted. 

The  aim  of  the  present  numerical  investigation  was  to  simulate  the  time  dependent 
liquid  flow  in  a  closed,  partially-filled  cylindrical  container.  The  spinrate  of  the  cylinder  was 
sufficiently  large  so  that  the  liquid  free  surface  was  also  cylindrical  and  centered  on  the  spin 
axis.  In  this  state  the  fluid  was  considered  to  be  in  solid  body  rotation.  The  spinrate  was 
then  perturbed  by  an  instantaneous  acceleration  to  a  new  steady  spin  resulting  in  a  time 
dependent  response  of  the  fluid  to  the  change  in  rotation  of  the  cylinder.  The  partially- 
filled  cylinder  problem  had  two  special  requirements.  The  computational  grid  had  to  be 
sufficient  to  resolve  a  possible  free  surface  shear  layer,  as  well  as,  boundary  layers  on  the 
solid  walls  of  the  container  (sidewall  and  endwalls).  In  addition,  the  boundary  conditions 
on  the  free  surface  had  to  be  posed  in  a  compatible  way  with  the  numerical  framework. 

In  this  study  the  Navier-Stokes  equations  were  written  in  cylindrical  coordinates. 
Since  the  symmetry  axis  of  the  cylindrical  container  was  identical  to  the  spin  axis  of  the  sys¬ 
tem,  the  formulation  was  written  assuming  an  axisymmetric  flowfield.  The  Navier-Stokes 
equations  were  written  in  stream  function- vorticity  form.  This  choice  permits  solution  of 
the  equations  using  standard  methods  without  the  use  of  artificial  compressibility. 

The  method  of  artificial  compressibility  was  developed  by  Chorin10  to  solve  the  prim¬ 
itive  variable  formulation  of  the  incompressible  Navier-Stokes  equations.  In  this  form  an 
inherent  numerical  difficulty  is  created  due  to  the  absence  of  a  time  derivative  in  the 
continuity  equation.  To  facilitate  the  construction  of  am  implicit,  approximately  factored 
numerical  method,  a  fictitious  time-dependent  term,  designed  to  vanish  as  steady  state  is 
approached,  is  used  to  modify  the  continuity  equation.  However,  approximately  factorized 
schemes  incur  errors  due  to  the  factorization.  When  this  type  of  scheme  is  combined  with 
the  artificial  compressibility  terms  in  the  continuity  equation,  the  momentum  equations 
are  then  contaminated  by  this  factorization  error. 

The  choice  of  nonprimitive  variables  allows  the  solution  of  the  governing  equations 
without  augmentation  with  artificial  terms.  The  Navier-Stokes  equations  were  discretized 
on  a  gild  utilizing  logarithmic  stretching  to  cluster  grid  points  near  solid  boundaries  and  the 
free  surface.  Second-order  finite-differences  were  used  for  all  derivative  terms  in  the  gov¬ 
erning  equations  except  the  convection  terms  for  which  upwind  differencing  vcas  employed. 
Finally,  the  equations  are  numerically  relaxed  using  a  Gauss-Seidel  iteration  method. 

II.  FREE  SURFACE  CONSTRAINTS 

For  this  numerical  investigation  the  angular  velocity  of  the  cylinder  is  considered  to 
be  large  enough  that  the  fluid  forms  an  annulus  of  near  uniform  thickness  on  the  vertical 
wall.  The  following  conditions  are  assumed: 

1.  Initially,  the  container  is  rotating  with  spin  rate  R  and  the  fluid  is  in  a  state  of  rigid- 

body  rotation  (i.e.  the  azimuthal  component  of  velocity  is  directly  proportional  to 

radial  position,  V  =  S7,/?  and  all  other  components  of  velocity  are  zero). 


2.  The  free  surface  is  cylindrical  and  centered  on  the  spin  axis.  This  implies  that  the 
radial  acceleration  is  much  greater  than  the  gravitational  acceleration  (5),  or  the 
rroude  number,  Fr  =  Q]a/g  1  for  a  cylinder  of  radius  a.  This  constraint  on  Q, 
results  in  a  Froude  number  that  is  larger  than  the  Ekman  number,  E  =  v/Qtc2  where 
c  is  the  cylinder  half-height  and  v  is  the  kinematic  viscosity  of  the  liquid. 


Winch11  and  Gerber12  derived  a  set  of  equations  that  describe  the  shape  of  the  free 
surface  in  a  partially-filled  cylinder  in  steady  state  axial  rotation.  These  equations  verify 
that  a  chosen  set  of  conditions  (Fr,  cylinder  aspect  ratio  and  fill  ratio)  will  result  in  a 
vertical  free  surface.  For  a  right  circular  cylinder  that  is  partially-filled  with  liquid  and 
that  is  suddenly  rotated  at  a  constant  angular  velocity  f 2/,  the  fluid  will  rise  along  the 
sidewall  of  the  cylinder  and  form  a  parabolic  free  surface.  When  the  fluid  finally  attains 
the  angular  velocity  of  the  cylinder,  a  steady  state  or  equilibrium  shape  is  reached  and  is 
defined  bv: 


Z  =  ZQ  + 


q}r 2 
~2g~ 


(1) 


Equation  1  describes  a  parabola  with  a  vertex  that  intersects  the  spin  axis  of  the 
cylinder  at  a  height  of  Z„.  Using  Equation  (1)  and  conservation  of  liquid  mass.  Gerber12 
has  derived  an  equation  for  the  shape  of  the  free  surface  when  the  free  surface  intersects 
both  endwalls  of  the  cylinder: 


Z  = 


Q)a2 
4  eg 


n2 


Z,  +  c+^(R 
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(2) 


where  c  is  the  cylinder  half-height,  a  is  the  radius,  and  Z,  is  the  initial  (£lj  =  0)  liquid 
level  in  the  cylinder. 


Xondimensionalizing  2  =  Z/a  and  r  =  R/a  and  defining  the  fill  ratio,  /  =  Z,/(2c), 
and  aspect  ratio,  a  =  c/a.  we  write  Equation  (2)  as, 


(3) 


When  the  free  surface  intersects  both  the  bottom  wall  (2  =  0)  and  the  top  wall 
2  =  2a)  of  the  cylinder,  we  can  define  the  radial  location  of  the  free  surface  at  2  =  0, 


r,  = 
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and  the  radial  location  of  the  free  surface  at  2  =  2a, 
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(5) 


For  a  given  set  of  parameters  (/  =  0.95.  a  =  3.0),  Figure  1  shows  the  difference  in 
free  surface  radial  location  between  the  cylinder  endwalls  (r2  —  rj)  as  a  function  of  Froude 
number.  The  free  surface  first  intersects  both  end  walls  at  about  Fr  =  109. 
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The  free  surface  is  essentially  vertical  for  Fr  above  900.  For  this  Fr  (a  =  0.1  meters), 
Q.)a  =  8S29  meters/sec2  >  g  =  9.81  meters/sec2  and  the  slope  of  the  free  surface  is 
nearly  infinite.  An  equation  for  the  radial  location  of  the  free  surface  can  be  derived  from 
conservation  of  liquid  mass13,  or  by  setting  Fr  =  oo  in  Equation  (4)  or  (5). 

r,  =(l-/),/2  (6) 

where  r,  =  0.2236  for  the  parameters  specified  above. 


III.  GOVERNING  EQUATIONS 

The  motion  of  fluid  within  a  rotating  cylinder  is  governed  by  the  unsteady,  incom¬ 
pressible  Navier-Stokes  equations.  The  natural  coordinate  system  for  this  problem  is  the 
cylindrical  system  (r,  #,  z)  where  u,r,and  id  are  the  radial,  azimuthal,  and  axial  veloc¬ 
ity  components  respectively.  Disturbances  to  the  motion  of  the  cylindrical  container  are 
axisymmetric,  thus,  the  flowfield  is  axisymmetric.  The  equations  are  written  in  axisvm- 
inetric  form  by  eliminating  all  8  derivatives.  The  dimensionless  form  of  the  equations  in 
the  nonrotating  reference  frame  are; 

Continuity  equation 


1  d(ru)  dw 
r  dr  dz 

Momentum  equations 


(7) 
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where, 

r  =  R/a,  z  -  Z/a ,  u  =  U/^lja,  v  =  V/Slja,  w  =  WjVlja 
p  =  P/p?l2a 2,  r  =  TQ/,  Re  =Q,ja2/v 

The  stream  function- vorticitv  form  of  the  Navier-Stokes  equations  is  used  in  the  so¬ 
lution  scheme.  The  definition  of  these  variables  are: 


Stream  function  ipz  =  ru 
rpr  =  —  rw 

Vorticity  C  =  uz  ~  wr 
Circulation  7  =  rv 


(11) 

(12) 

(13) 

(14) 


The  governing  equations  cast  in  this  form  are  obtained  as  follows.  The  first  equation 
of  the  new  system  is  derived  from  the  definition  of  vorticity  (Eq.  13)  using  the  definition 
of  stream  function  (Eqs.  11,12).  This  equation  is  elliptic  and  is  in  the  form  of  a  Poisson 
equation.  The  radial  momentum  equation  (Eq.  8)  is  differentiated  with  respect  to  z 
and  from  this  equation  the  r  derivative  of  the  axial  momentum  equation  (Eq.  10)  is 
subtracted.  This  cross-differentiation  eliminates  pressure  and  yields  a  parabolic  vorticity 
transport  equation.  The  azimuthal  momentum  equation  (Eq.  9)  is  recast  into  a  parabolic 
equation  using  the  definition  of  circulation  (Eq.  14).  The  stream  function-vorticity  form 
of  the  Navier-Stokes  equations  are: 


d2i'  d24 ' 


dr2 


+ 


1  drl> 

-;a7  =  rC 


<3C  ,  ,  <9C 

dr  +  dr  +  Dz 


u  27  dy  _  1 
r  r3  dz  Re 


n  5  !?( 

dr 2  dz2  r  dr 


(15) 

(16) 


and  the  circulation  equation  is: 


dl  ,  d 7 

dr  dr  dz 


d2-y  d2  7 

d^2  +  Ik2 


(17) 


Equation  16  is  nonlinear  in  the  convective  terms  (§£  +  wjfc'j  since  u  and  w  are  functions  of 
the  dependent  variable  (  via  Equation  13.  The  equation  is  parabolic  in  time  and  poses  an 
initial  value  problem.  The  same  is  true  of  the  circulation  equation  (Eq.  17).  The  stream 
function  equation  (Eq.  15)  is  elliptic  and  poses  a  boundary  value  problem  which  can  be 
solved  by  iterative  methods.  Equation  15  is  de-coupled  from  the  vorticity  equation  (Eq. 
16)  and  the  circulation  equation  (Eq.  17).  As  a  result,  Equations  16  and  17  can  be  solved 
first,  followed  by  the  separate  solution  of  Equation  15  with  vorticity  as  a  source  term. 

IV.  BOUNDARY  CONDITIONS 

The  governing  equations  are  applied  to  the  domain  of  interest  as  shown  in  Figure 
2.  In  the  case  of  a  cylinder  full  of  liquid  the  boundaries  of  the  domain  are  the  endwalls 
(Z  =  0  and  Z  =  2c),  the  sidewall  (R  =  a)  and  the  geometric  axis  ( R  =  0)  about  which 
the  cylinder  spins.  For  the  partially-filled  cylinder,  the  boundary  at  R  —  0  is  replaced  by 
a  free  surface  located  at  R  —  R,.  The  location  of  the  free  surface  was  discussed  in  Section 
II  and  given  by  Equation  6.  The  boundary  conditions  for  the  free  surface  will  be  discussed 
later  in  this  section. 
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The  governing  equations  are  written  in  the  inertial  (nonrotating)  reference  frame,  thus, 
the  no-slip  boundary  conditions  on  the  cylinder  walls  are:  U  =  W  =  0,  V  =  Slja  for  R  =  a; 

U  =  W  =  0,  V  =  QjR  for  Z  =  0;  and  U  ~  W  =  0,  V  =  ft/R  for  Z  =  2c.  These  no-slip 

boundary  conditions  are  made  dimensionless  using  the  same  definitions  employed  for  the 
governing  equations  and  become: 

r  =  1  :  u  =  =  0,  t>  =  1  (18) 

z  =  0'-u  =  w  =  0,v  =  r  (19) 

z  =  2a  :  u  =  w  =  0,v  =  r  (20) 

The  boundary  conditions  for  the  cylinder  axis  are  derived  from  the  nondimensional 
continuity  equation  in  primitive  variable  form  (Eq.  7): 

1  d(ru)  dw  n  du  dw  .  . 

-r—  +  S7=0  or  ra7  +  “  +  raT  =  0  (21) 

Now,  setting  r  =  0,  we  obtain  u  =  0.  Multiplying  both  sides  of  the  axial  momentum 
equation  (Eq.  10)  by  r  and  setting  r  =  0  we  obtain  the  condition  dw/dr  =  0.  Multiplying 
both  sides  of  the  azimuthal  momentum  equation  (Eq.  9)  by  r2  and  setting  r  =  0we  obtain 
the  condition  v  =  0.  Hence, 


r  =  0:u  =  u  =  0,  ——  =  0  (22) 

dr 

Because  of  flow  symmetry  we  restrict  the  domain  bounded  by0<r<l,0<2<2a 
to  a  computational  domain  bounded  by0<r<l,0<z<a  (Figure  2).  In  Reference  13, 
computational  verification  of  flow  symmetry  was  performed  over  the  domain  2  =  0  to  2a 
for  a  cylinder  full  of  liquid.  We  term  this  situation  the  “gyroscope  problem”  since  the 
entire  cylinder  undergoes  the  same  motion.  Once  the  flowfield  for  0<r<l,0<z<a 
(one  quadrant)  is  computed,  the  flowfield  in  the  rest  of  the  cylinder  can  be  obtained. 
This  domain  restriction  creates  an  additional  boundary,  the  symmetry  plane,  requiring  a 
boundary  condition.  On  the  symmetry  plane  we  require  no  axial  flow  and  no  axial  fle  w 
gradients  for  the  radial  and  azimuthal  velocity  components.  Thus, 

du  dv 

z  =  a:u;  =  0,^  =  0,—  =  0  (23) 

The  above  restriction  of  the  domain  using  symmetry  arguments  does  not  apply  when 
different  parts  of  the  cylinder  have  different  angular  velocities  as  in  the  “centrifuge  prob¬ 
lem”.  In  this  case,  the  bottom  wall  and  the  side  wall  of  the  cylinder  have  an  angular 
velocity,  ft/,  that  is  not  the  same  as  the  angular  velocity  of  the  top  wall,  ft^.  As  a  result 
two  quadrants  (0<r<l,0<z<  2a)  of  the  flowfield  must  be  considered  as  the  flow 
is  symmetric  about  the  spin  axis  but  not  the  midplane  of  the  cylinder.  This  adds  the 
boundary  condition: 

ft', 

z  =  2a:i£  =  u;  =  0,t;  =  zzr~r  (24) 

S  Ij 

The  boundary  conditions  are  written  in  the  variables  of  stream  function,  vorticity, 
and  circulation.  For  the  gyroscope  problem  the  conversion  of  boundary  conditions  to  the 
nonprimitive  variables  is  as  follows: 


mmsmm 


Axis: 

r  =  0 

with  Equations  11  and  12  yields: 

*0  =  0 

(r  =  0) 

U  =  IL'r  = 

=  0  with  Equation  13  yields: 

C  =  0 

v  —  0 

with  Equation  14  yields: 

7  =  0 

*P  -  0,  C  =  0,  7  =  0 

(25) 

Sidewall: 

u  =  w  = 

0  with  Equations  11  and  12  yields: 

0  =  0 

(r  =  1) 

u2  =  0 

with  Equation  13  yields: 

C  =  ~Wr  =  VVr 

V  =  1 

with  Equation  14  yields: 

7  =  1 

4'  =  0,  (  =  4>rr,  7=1 

(26) 

End  wall: 

u  —  w  = 

0  with  Equations  11  and  12  yields: 

0  =  0 

(*  =  o) 

wr  =  0 

with  Equation  13  yields: 

C  =  U2  =  7022 

v  =  r 

with  Equation  14  yields: 

7  =  r2 

4>  =  0,  C  =  -4’ 22,  7  =  r2 

r 

(27) 

Svmmetrv  Plane: 
(z  =  a ) 

w  =  0 

with  Equations  11  and  12  yields: 

0  =  0 

Uj  =  w  = 

=  0  with  Equation  13  yields: 

C  =  0 

r2  =  0 

with  Equation  14  yield?: 

0 

11 

N 

0  =  0,  C  =  0,  72  =  0 


(28) 


For  the  centrifuge  problem  the  fluid  is  driven  by  the  top  wall  of  the  cylindrical  con¬ 
tainer.  This  top  wall  has  a  steady  spinrate  that  is  not  the  same  as  the  spinrate  for  the  rest 
of  the  cylinder.  As  noted  in  the  derivation  of  the  governing  equations,  the  spinrate  used 
to  yield  dimensionless  equations  is  Slj.  The  boundary  condition  on  the  symmetry  plane 
(Eq.  28)  is  replaced  by: 


End  wall  (z  =  2a)  :  0  =  0, 


C  =  -0, 

r 


2 

1  =  W 


(29) 


The  free  surface  boundary  conditions  on  7  and  £  can  be  derived  from  the  assumption 
that  the  free  surface  does  not  support  tangential  stress  (i.e.  r$r  =  tzt  =  0). 


h 

d(v/r)  1 

r  _  4- 

_  h 

d(vr/r 2) 
r - 

_  ^ 

'a(7/r2)' 

2 

3r  r  5# 

2 

3r 

2 

dr 

=  0 


or 


d(7/r2) 

dr 


=  0  (30) 
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T‘t  2 


3u>  P  du  dw  5u  u  f  ,  2  ,  1  n  _  ,  , 

+  ~k~  =  ~  —  ^  h  I-  2—  =  —  — C  H  ^2  —  0  or  2?/’22  —  r(  —  0 
dr  J  2  02  dr  d2  2  1  r 


Since  the  stream-function  equation  (Eq.  15)  is  elliptic  and  its  solution  represents  a  bound¬ 
ary  value  problem,  the  boundary  condition  for  stream  function  on  the  free  surface  should 
be  a  Dirichlet  or  Neumann  condition.  Since  the  condition  ip  =  0  is  imposed  on  the  other 
boundaries  (endwalls  and  sidewall)  of  the  cylinder,  the  Dirichlet  condition  ip  =  0  is  also 
imposed  on  the  free  surface.  This  condition  is  equivalent  to  setting  u  —  0  on  the  free 
surface  and  is  consistent  with  the  vertical  free  surface  assumption  made  in  Section  II.  As 
a  result,  the  free  surface  is  characterized  as  a  “slippery  rod”  of  radius  r,.  The  free  surface 
(r  =  r,)  boundary  conditions  are  prescribed  as: 

*  =  °.  <  =  0-  f(^)  =  0  (32) 


From  Equation  32,  values  of  7  on  the  free  surface  can  be  obtained  from  values  of  7  in 
the  interior  flow  using  a  second  order  forward  finite- difference  in  the  radial  direction.  Use 
of  Equation  32  as  boundary  conditions  for  the  free  surface  is  restrictive.  Surface  waves  and 
free  surface  shape  changes  are  ignored.  Unfortunately,  these  phenomena  have  been  shown 
to  exist  for  rotating  and  nutating  free  surface  flows.  Stewartson1,  in  predicting  the  flow  in 
a  partially-filled  cylinder,  derived  a  free  surface  boundary  condition  which  included  small 
perturbations  to  the  free  surface  shape  (r,  +  v)-  In  the  present  formulation,  Stewartson’s 
condition  yields  a  time-dependent  boundary  condition  for  ip  on  the  free  surface.  However, 
since  Equation  15  consists  of  only  spatial  derivatives  in  ip,  the  boundary  condition  on  ip 
cannot  be  time-dependent.  Therefore,  Stewartson’s  condition  cannot  be  used  with  the 
incompressible  Navier-Stokes  equations  in  stream  function-vorticity  form. 


Free  surface  boundary  conditions  such  as  Equation  32  were  used  by  Harlow  and 
Welch14  in  a  paper  that  described  the  marker  and  cell  method.  In  this  work  the  marker  and 
cell  technique  was  applied  to  the  time-dependent  flow  of  an  incompressible  fluid,  the  bound¬ 
ary  of  which  was  partially  confined  and  partially  free.  The  time-dependent  Navier-Stokes 
equations  were  written  in  primitive  variable  form  and  free  surface  boundary  conditions 
were  required  for  all  variables.  Velocity  boundary  conditions  at  a  free  surface  were  derived 
from  the  continuity  equation  while  the  pressure  boundary  condition  was  based  on  a  con¬ 
stant  external  (or  applied)  pressure.  They  observed  that  the  pressure  condition  should  be 
based  on  the  requirement  of  vanishing  normal  stress  component  at  the  free  surface,  but 
this  would  be  difficult  to  implement  unless  the  orientation  of  the  surface  was  known  at  all 
times.  They  observed  that  such  a  determination  was  unreasonably  difficult  for  a  finite- 
difference  scheme  (assuming  an  Eulerian  approach)  and  equated  the  pressure  on  the  free 
surface  to  the  applied  pressure.  Harlow  and  Welch  achieved  remarkably  accurate  results 
for  a  wide  class  of  free  surface  problems  in  uncontained,  nonrotating  flow  systems. 


V.  NUMERICAL  SOLUTION  SCHEME 


This  section  gives  the  numerical  solution  scheme  for  the  governing  equations  with 
the  boundary  conditions  presented  in  preceding  sections.  First  an  overview  of  the  solution 


methodology  is  presented  which  is  similar  to  those  used  in  previous  references.15'16'17’18  The 
computational  domain  for  the  gyroscope  problem  is  bounded  by  one  cylinder  endwall,  the 
symmetry  plane,  the  cylinder  sidewall  and  the  free  surface  (Figure  2).  For  the  centrifuge 
problem  the  domain  is  bounded  by  the  cylinder  endwalls,  the  cylinder  sidewall  and  the 
free  surface.  A  grid  is  constructed  in  the  domain,  and  the  solution  is  computed  at  mesh 
points  on  the  grid.  The  computation  begins  with  initial  values  for  t/>,  (,  7,  u,  u,  and  w  at 
every  mesh  point  and  at  time  equal  to  zero.  For  a  time  accurate  solution  where  transient 
flow  phenomenon  are  of  interest,  the  initial  solution  corresponds  to  a  real  flow  condition: 

=  0,  C  =  0,  7  =  77- r2,  u  =  0,  v  =  ^-r,  w  =  0  (33) 

i  1/ 

At  time  level  r  =  r  +  Ar,  the  computational  cycle  begins  with  the  simultaneous  solution  of 
the  finite-difference  analog  of  Equations  16  and  17  for  £  and  7  at  all  interior  mesh  points. 
The  governing  equations  are  discretized  directly  on  a  stretched  grid.  This  procedure  has 
been  used  by  other  researchers6.  The  next  step  is  to  solve  the  finite  difference  analog  of 
Equation  15  at  all  interior  mesh  points.  This  equation  involves  computed  values  of  £  as  a 
source  term.  These  equations  are  solved  by  Gauss-Seidel  relaxation  using  successive  over¬ 
relaxation  (SOR).  At  this  point  new  velocity  components  are  evaluated  from  the  definitions 
of  stream  function  (Equations  11  and  12)  and  circulation  (Equation  13).  Finally,  the 
boundary  conditions  are  imposed.  Iterative  solution  of  the  flowfield  for  the  current  time 
step  is  continued  until  convergence.  The  computational  cycle  is  repeated  until  the  desired 
time  is  reached. 


£ 

& 


Stretching  of  the  grid  in  the  computational  domain  is  accomplished  with  the  following 
equation  (after  Roberts19): 


r  = 


md  +  l)/(d  -  1)]"  -  1} 
1  +  l(d+l)/(d-  !)]'> 


(34) 


where  a  uniformly  spaced  coordinate,  /?(x),  can  be  transformed  into  a  nonuniform  coor¬ 
dinate,  r(x),  that  varies  continuously  throughout  the  domain  but  is  clustered  near  the 
boundaries.19  The  parameter  d  <  1  is  the  control  to  achieve  desired  clustering  with  very 
fine  boundary  spacing  for  d  <  1,  Similarly  for  the  axial  direction,  a  uniformly  spaced 
coordinate,  tj(x ),  can  be  transformed  into  a  nonuniform  coordinate,  z(x),  using: 


z  = 


(35) 


q(l  -6  +  (l  +  6)[(6  +  l)/(6  -  l)]"-1) 

1  +  K& +!)/(* -I)]-’1 

The  parameter  b  <  1  is  defined  in  the  same  manner  as  d.  These  formulas  are  derived  in 
Reference  13.  For  a  full  cylinder,  grid  points  are  clustered  at  R  =  a  and  Z  =  0.  For  a 
partially-filled  cylinder,  grid  points  are  clustered  at  R  =  a,  R  =  Rt  and  Z  =  0.  Figures  3 
and  4  show  typical  grids  for  the  full  and  partially-filled  cases  ( c/a  =  1.0). 

The  finite-difference  analog  of  the  governing  equations13  are  obtained  using  central  and 
upwind  difference  expressions  for  the  derivative  terms.  The  convective  terms  of  Equations 
16  and  17  are  represented  using  upwind  differencing  (see  Appendix  C,  Reference  13).  All 
other  derivative  terms  are  represented  using  central  differences.  All  spatial  derivatives  are 
represented  by  three-point  second  order  finite-difference  formula.  The  time  derivative  is 
represented  by  a  second  order  difference  formula. 
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The  iterative  solution  of  the  governing  equations  is  accomplished  using  Gauss-Seidel 
relaxation.17’20  Gauss-Seidel  iterations  converge  very  rapidly  provided  that  a  condition 
of  diagonal  dominance  exists  for  the  coefficient  matrix  of  the  system  of  equations.  For 
the  governing  equations  formulated  in  stream  function,  vorticity  and  circulation  variables, 
diagonal  dominance  has  been  shown  to  exist.15  The  present  study  uses  a  modification  of  the 
Gauss-Seidel  method  called  successive  over-relaxation  (SOR).  This  technique  accelerates 
the  convergence  of  the  algorithm  and  as  a  result  reduces  the  required  computational  time. 
As  the  Gauss-Seidel  iteration  is  applied  to  the  system  of  simultaneous  algebraic  equations 
several  iterations  are  made  before  convergence  to  an  acceptable  level  is  achieved.  During 
this  process  we  observe  the  direction  of  the  change  in  the  unknown  at  a  grid  point  between 
two  successive  iterations  and  anticipate  that  the  same  trend  will  continue  on  to  the  next 
iteration.  To  accelerate  convergence,  an  arbitrary  correction  to  the  intermediate  values  of 
the  unknown  is  made: 

Ui?1'  =  «£  +  w  (uf,+1  -  «*,')  (36) 

Here,  k  denotes  iteration  level  and  u*+1  is  the  most  recent  value  of  calculated  from 
the  Gauss-Seidel  procedure,  is  the  value  from  the  previous  iteration  as  adjusted  by 
previous  application  of  this  formula,  and  uf*1’  is  the  guess  for  Uij  at  the  k  +  1  iteration 
level.  The  formula  is  applied  immediately  at  each  point  after  u^+1  has  been  obtained  and 
uf/1  replaces  u*+1  in  all  subsequent  calculations  in  the  cycle.  The  relaxation  parameter  is 
u>.  For  over- relaxation  1  <  u  <  2.  For  convergence  u>  <  2.  In  the  present  study  a  constant 
value  of  1.75  was  used. 


VI.  COMPARISONS  WITH  PREVIOUS  WORK 


The  computational  scheme  is  initially  checked  for  accuracy  by  comparison  with  flow 
velocity  data  from  experiments  with  a  filled  right  circular  cylinder.  Watkins  and  Hussey21 
used  LDV  to  measure  the  azimuthal  velocity  component  in  a  cylinder  of  aspect  ratio  unity, 
spinning  up  to  an  angular  rate  of  1.83  radians  per  second  from  rest.  The  Reynolds  number 
of  the  flow  was  9741.6  based  on  the  final  spinrate  of  the  container.  Azimuthal  velocity 
was  measured  radially  across  the  midplane  of  the  cylinder  at  various  times  in  the  spin- 
up  process.  This  is  compared  to  the  present  work  in  which  a  31  x  31  (r,  z)  grid  was 
used  with  logarithmic  stretching  in  the  radial  ( d  =  0.05)  and  axial  (6  =  0.05)  directions. 
Figure  3  illustrates  the  computational  grid.  The  solution  was  started  with  the  cylinder 
at  rest  and  instantaneously  given  the  spinrate  of  1.83  rad/sec.  This  process  approximates 
the  experimental  procedure  in  which  the  container  required  a  small  amount  of  time  to 
achieve  full  rotational  speed.  Watkins  and  Hussey  recorded  velocity  profile  measurements 
at  various  time  intervals  between  42.8  and  210  seconds  during  the  spin-up  process.  The 
computation  used  a  small  numerical  time  step  (0.0183  seconds)  to  insure  numerical  stability 
and  was  checked  at  45,  75,  105  and  150  seconds  in  the  spin-up  calculation.  Figure  5  shows 
a  comparison  between  the  measured  and  computed  azimuthal  velocity  profile  across  the 
midplane  of  the  cylinder  at  these  times.  In  addition,  the  solution  by  Kitchens4  is  shown. 
In  all  cases  the  comparison  to  measured  azimuthal  velocity  is  good.  Completion  of  the 
spin-up  process  is  observed  after  150  seconds  as  the  azimuthal  velocity  profile  takes  on  a 
linear  behavior  (V  =  Cl/R). 
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A  second  experimental  case  by  Watkins  and  Hussey  provides  an  additional  comparison 
between  the  numerical  code  and  LDV  measurements  for  a  cylinder  full  of  liquid.  For  this 
case  the  aspect  ratio  of  the  cylinder  was  increased  to  1.515  and  the  Reynolds  number 
was  decreased  to  3076.  The  identical  grid  and  logarithmic  stretching  constants  were  used. 
Again  the  cylinder  was  initially  at  rest  and  was  instantly  given  an  axial  spinrate  of  1.34 
radians  per  second.  The  computational  time  step  was  again  chosen  as  0.0183  seconds 
and  the  computation  was  checked  at  15,  45,  75,  and  120  seconds.  Figure  6  shows  good 
comparisons  between  measured  and  computed  azimuthal  velocity  profiles  on  the  midplane 
of  the  cylinder.  For  this  case  the  smaller  Reynolds  number  indicates  a  shorter  spin-up  time 
and  this  is  indeed  observed. 

A  third  computation  is  compared  to  the  solution  of  the  von  Karman  problem.  This 
is  the  steady  flow  of  an  incompressible  viscous  liquid  due  to  an  infinite  rotating  disk.  A 
stationary  liquid  occupies  the  semi-infinite  region  on  one  side  of  the  disk,  the  motion  of 
which  is  rotationally  symmetric.  The  effect  of  the  disk  is  to  move  the  fluid  near  its  surface 
radially  outwards,  and  in  turn  induce  an  axial  flow.  The  main  interest  in  this  problem  is 
that,  by  virtue  of  assumptions  about  the  velocity  components,  the  Navier-Stokes  equations 
reduce  to  a  set  of  ordinary,  nonlinear  differential  equations  in  a  single  independent  variable. 
Several  papers  have  been  devoted  to  the  solution  of  this  set  of  equations,  including  Rogers 
and  Lance22. 

For  a  cylinder  of  aspect  ratio  unity,  initially  at  rest  and  instantaneously  accelerated 
to  an  axial  spinrate  of  1.83  radians  per  second,  the  resulting  Reynolds  number  is  9741.6. 
Early  in  the  spin-up  process  the  core  flow  is  not  rotating  (azimuthal  velocity  is  zero),  and 
only  the  fluid  near  the  endwalls  and  sidewall  is  effected  by  the  container’s  motion.  The 
present  computation  was  performed  using  a  time  step  of  0.0055  seconds  and  marched  to  a 
time  of  4.54  seconds.  Figure  7  shows  the  azimuthal  velocity  profile  very  near  the  endwall 
( Z/a  =  0.092)  and  approximately  outside  the  Ekman  layer  (for  this  case  the  Ekman  layer 
thickness  is  of  order  0.01  or  E1^2).  Note  that  the  flow  between  the  cylinder  axis  and  about 
R/a  =  0.82  is  not  rotating.  This  region  of  fluid  meets  the  criteria  for  the  von  Karman 
problem,  i.e.,  a  rotating  disk  beneath  a  fluid  at  rest.  Figure  8  shows  the  radial  velocity 
profile  at  a  radius  of  R/a  =  0.76,  near  the  cylinder  endwall  and  in  the  Ekman  layer.  The 
agreement  with  Rogers  and  Lance  and  the  Navier-Stokes  solution  of  Kitchens4  is  good. 
The  radial  velocity  profile  crosses  into  values  less  than  zero  at  about  Z/a  —  0.07.  This 
is  generally  defined  as  the  edge  of  the  Ekman  layer  because  the  radial  velocity  profile 
passes  through  zero  at  one  axial  point  for  a  particular  radial  location  (exceptions  occur  for 
transient  reverse  flow  regions).  The  initial  guess  for  the  order  of  thickness  for  the  Ekman 
layer  (0.01)  is  also  verified. 

Another  check  on  the  computation  comes  from  inertial  oscillations  which  occur  when 
the  angular  velocity  of  the  cylinder  is  suddenly  increased.  These  inertial  oscillations  are 
transients  in  the  flowfield  that  eventually  cease.  Warn-Varnas  et  al23  used  LDV  to  detect 
the  presence  of  these  oscillations  in  a  cylinder  of  aspect  ratio  0.3182,  full  of  water  and 
disturbed  from  an  initial  state  of  solid  body  rotation  at  0.6282  radians  per  second  to  a 
new  spinrate  of  0.7678  radians  per  second.  The  resulting  Reynolds  number  of  the  flow 
(based  on  the  final  spin  rate)  was  7334.  They  collected  azimuthal  velocity  data  at  several 
locations  in  the  flowfield  during  the  spin-up  of  the  liquid.  These  velocity  measurements 


were  then  reduced  to  a  nondimensional  zonal  velocity  defined  by, 

Zonal  Velocity  =  £2/(  1  —  v/r)/(f2/  —  £2,)  (37) 

which  takes  on  the  value  of  zero  when  the  fluid  has  attained  a  new  state  of  solid  body 
rotation  (v  =  r).  The  present  computation  was  performed  for  this  case  using  a  43x21 
stretched  grid  and  a  time  step  of  0.0052  seconds.  Azimuthal  velocity  was  collected  for 
every  time  step  at  4  grid  points,  corresponding  to  the  locations  used  by  Warn-Varnas. 

Figures  9  thru  12  show  the  comparisons  between  measured  and  computed  zonal  veloc¬ 
ity  versus  time,  for  two  axial  locations  at  a  radial  location  of  R/a  =  0.25  (in  the  core  flow 
near  the  spin  axis  of  the  cylinder).  Figure  9  shows  an  axial  location  of  Z/a  =  0.3182  (at 
the  midplane  of  the  cylinder)  and  Figure  10  shows  an  axial  location  of  Z/a  =  0.079  (near 
the  cylinder  endwall).  The  Navier-Stokes  solution  of  Kitchens4  is  also  shown  in  Figure  9. 
Correlations  between  computation  and  measurement  are  within  the  experimental  accuracy 
reported  by  Warn-Varnas.  These  plots  illustrate  the  axial  structure  of  the  inertial  oscilla¬ 
tions.  The  amplitude  of  the  oscillations  are  nearly  the  same  at  both  axial  locations,  but  a 
phase  shift  has  occurred.  Figures  11  and  12  show  the  comparisons  between  measured  and 
computed  zonal  velocity  versus  time  for  the  same  axial  locations,  but  at  a  radial  location 
of  R/a  =  0.75  (near  the  cylinder  sidewall).  Once  again  the  axial  structure  of  the  inertial 
oscillations  are  observed;  however,  the  amplitude  of  the  oscillations  is  severely  damped. 
This  more  gradual  transition  to  solid  body  rotation  may  be  due  in  part  to  the  proximity 
of  these  points  to  the  sidewall  boundary  layers. 

In  the  case  of  a  partially-filled  cylinder,  the  free  surface  boundary  conditions  can  be 
checked  against  available  data  from  Shadday.9  He  measured  the  axial  and  azimuthal  veloc¬ 
ity  components  in  a  partially-filled  cylinder  initially  in  a  state  of  solid  body  rotation  and 
disturbed  by  the  acceleration  of  the  top  'wall  to  a  new  steady  state  spinrate.  The  rotation 
rate  of  the  cylinder  sidewall  and  bottom  wall  was  maintained.  For  this  configuration,  the 
flow  is  no  longer  symmetric  about  the  cylinder  midplane;  thus,  the  entire  flowfield,  for  one 
azimuthal  plane,  must  be  computed. 

The  cylinder  was  36.7%  full  of  liquid  ten  times  more  viscous  than  water.  In  the  exper¬ 
iment  the  entire  cylinder  was  rotated  at  1000  rpm  until  a  steady  state  flow  condition  was 
reached.  Then  the  top  wall  of  the  cylinder  was  accelerated  by  5%  in  spinrate  while  the  rest 
of  the  cylinder  was  maintained  at  1000  rpm.  The  system  was  maintained  until  steady  flow 
existed.  Based  on  a  spinrate  of  1000  rpm,  the  Reynolds  number  was  93518.0  (the  Ekman 
number  was  1.0119xl0“5).  The  aspect  ratio  of  the  cylinder  was  1.028.  LDV  measurements 
of  flow  velocity  profiles  were  made  radially  across  the  flowfield  at  four  specified  heights 
from  the  bottom  wall  (Figure  13). 

For  the  computation  of  this  flowfield,  the  free  surface  was  assumed  to  be  vertical 
and  located  at  a  nondimensional  radius  of  0.7957.  This  radial  location  was  measured  by 
Shadday  continuously  during  the  test  and  reported  to  be  constant.  As  a  verification,  the 
method  described  in  Section  II  was  used  to  check  the  shape  of  the  free  surface  under  these 
conditions.  For  the  axial  spinrate  in  this  experiment,  the  slope  of  the  free  surface  was 
0.690  degrees. 

The  grid  chosen  for  the  computation  consisted  of  30  points  m  *he  radial  direction  and 
61  points  in  the  axial  direction  (Figure  14).  Grid  stretching  was  used  to  cluster  grid  points 


near  the  solid  boundaries  and  the  free  surface.  The  grid  spacing  at  the  cylinder  sidewall 
and  endwalls  was  .0014  and  .0018,  respectively.  This  results  in  at  least  2  grid  points  within 
the  wall  boundary  layers.  The  thickness  of  the  sidewall  and  endwall  boundary  layers  are 
0(E V4)  and  0(E V2),  respectively.  Grid  clustering  on  the  free  surface  was  the  same  as 
that  used  on  the  cylinder  sidewall.  Several  grid  sensitivity  studies  were  conducted  to  insure 
that  the  chosen  grid  was  adequate. 

The  LDV  measurements  of  radial  velocity  profiles  were  recorded  by  Shadday  after 
the  top  wall  was  accelerated  to  a  new  spinrate  and  the  resultant  flow  had  achieved  steady 
state.  Shadday  reports  a  predicted  time  of  6.0  seconds  for  the  flow  to  reach  steady  state 
((c2 Shadday’s  time- dependent  computation  was  run  for  7.4  seconds.  The  cur¬ 
rent  time- dependent  solution  was  also  run  for  7.4  seconds  and  studies  were  made  to  verify 
that  the  solution  had  indeed  reached  a  steady  state.  It  was  observed,  however,  that  steady 
state  was  reached  well  before  6.0  seconds. 

Figures  15  through  22  show  the  comparison  between  measured  and  computed  velocity 
components  at  axial  stations  of  Z/a  =  0.497,  0.9,  1.304,  and  1.706.  In  each  case,  the  nondi- 
mensional  axial  and  azimuthal  velocity  is  hown  as  a  function  of  radial  position.  Figure  13 
illustrates  the  locations  of  these  axial  stations  relative  to  the  differentially  rotating  end- 
wall.  For  the  azimuthal  velocity,  the  agreement  between  measurements  and  computations, 
using  both  the  present  code  and  Shadday’s  results,  is  good  for  all  axial  stations.  For  the 
axial  velocity,  agreement  with  measurements  for  both  codes  is  good  except  near  the  free 
surface.  The  magnitude  of  this  discrepancy  is  nearly  constant  for  all  axial  stations. 

The  discrepancy  between  computed  and  measured  axial  velocity  near  the  free  surface 
may  be  due  in  part  to  a  velocity  bias  in  the  LDV  system,  as  reported  by  Shadday.  This 
discrepancy  may  also  be  due  to  the  free  surface  shape  and  boundary  condition  assumptions. 
Recall  that  the  present  work  assumes  that  the  free  surface  is  vertical.  An  attempt  was 
made  to  minimize  this  effect  in  the  experiment  by  choosing  an  axial  spinrate  that  yielded  a 
free  surface  slant  of  less  that  a  degree.  In  addition,  the  boundary  conditions  were  obtained 
using  the  assumption  of  inviscid  free  surface  flow  (zero  tangential  shear  stress)  and  constant 
normal  pressure  on  the  free  surface.  Although  the  effect  of  the  inviscid  boundary  condition 
on  the  final  result  is  not  known,  this  boundary  condition  would  be  violated  most  severely 
near  the  differentially  rotating  endwall. 

VII.  RESULTS  AND  DISCUSSION 

In  this  section,  the  resultant  roll  moment  exerted  on  the  cylinder  is  examined  for  full 
and  partially-filled  cylinders  in  rapid  axial  rotation.  The  container  is  initially  in  a  state  of 
solid  body  rotation  with  an  axial  spin  rate  that  is  suddenly  accelerated.  The  parameters 
chosen  for  this  study  were  Reynolds  number  1000,  cylinder  aspect  ratio  3.0,  initial  spinrate 
300.0  rad/sec,  final  spinrate  400.0  rad/sec,  cylinder  radius  0.1  meters,  liquid-fill  density 
1400  kilograms /meter3  and  liquid  fill  ratios  of  .95  and  1.0. 

For  both  the  full  and  partially-filled  cases,  a  21  x  41  grid  was  selected  with  the  loga¬ 
rithmic  grid  clustering  at  the  endwalls,  sidewall  and  free  surface.  The  Ekman  number  for 


these  flow  conditions  is  l.llxlO-4  and  provides  guidance  for  the  selection  of  grid  clustering 
at  the  How  boundaries.  For  the  sidewall  and  free  surface,  the  grid  spacing  is  selected  as 
.0042  (<  E1^).  For  the  endwall,  the  grid  spacing  is  selected  as  .0089  (<  E 1/2).  Figures 
23  and  24  display  the  computational  grids  for  cylinders  with  fill  ratios  of  1.0  and  0.95, 
respectively.  The  location  of  the  free  surface,  r,  =  0.2236,  was  found  from  Equation  6. 
Figure  1  shows  the  difference  in  free  surface  radial  location  at  z  =  6  and  2  =  0,  A r,  as  a 
function  of  Froude  number.  For  the  initial  spinrate  of  300  rad/sec  (Fr  =  917.43,  a  =  0.1 
meters),  the  slope  of  the  free  surface  is  about  0.20  degrees.  The  numerical  computation 
was  carried  out  for  1000  time  steps  or  to  0.0025  seconds.  At  this  time,  the  liquid  has  suf¬ 
ficiently  adjusted  to  the  cylinder  spinrate;  and  the  resultant  liquid  roll  moment  has  been 
significantly  reduced. 


An  important  calculation  for  this  problem  is  the  roll  moment  exerted  on  the  cylinder 
by  the  liquid  fill.  The  moment  is  obtained  by  integrating  the  shear  stress  tangent  to 
the  cylinder  endwalls  and  sidewall  over  the  entire  cylinder.  The  contribution  to  the  roll 
moment  at  each  grid  point  on  the  cylinder  endwalls  is  proportional  to  rz6  which  for  the 
case  of  axisymmetric  flow  is  n(dv/dz).  The  contribution  to  the  roll  moment  at  each  grid 
point  on  the  cylinder  sidewall  is  proportional  to  rTg  which  is  equal  to  fj(dv/dr  —  v/r).  The 
roll  moment  is  made  dimensionless  using  a  convention  suggested  by  Murphy24: 


Clrm  = 


Roll  Moment 

(liquid  mass)( cylinder  radius)2 (final  spinrate)2 


where  the  liquid  mass  is  defined  by  the  container  volume  and  liquid  density.  The  liquid  roll 
moment  is  less  than  zero  because  the  liquid  acts  to  despin  the  container.  When  the  liquid 
is  in  a  state  of  solid  body  rotation,  before  the  perturbation  to  axial  spinrate  and  after  the 
liquid  has  adjusted  to  the  perturbation,  Cirm  is  zero.  This  can  be  used  as  an  indication 
of  the  liquid  spin-up  time.  Figure  25  shows  the  time  history  of  the  total  coefficient  of 
liquid  roll  moment.  Clrm  for  the  partially-filled  case  is  less  than  that  for  the  full  case. 
After  a  very  short  time,  Clrm  approaches  zero.  The  smaller  Clrm  for  the  partially-filled 
cylinder  suggests  that  either  the  wall  shear  stress  is  different  in  the  partially-filled  case  or 
that  the  part  of  the  cylinder  endwall  between  r  =  0  and  r  =  r,  for  the  full  case  is  a  major 
contributor  to  the  total  moment. 


Figures  26  and  27  show  the  sidewall  and  endwall  components  of  Clrm ■  The  sidewall 
contributes  significantly  to  the  total  roll  moment  and  shows  a  marked  difference  for  full 
and  partially-filled  cylinders.  For  the  cylinder  endwall,  difference  in  Clrm  for  the  partially- 
filled  and  full  cases  is  insignificant.  Figure  28  shows  the  azimuthal  velocity  profiles  near  the 
cylinder  sidewall  for  two  axial  locations.  These  values  are  computed  for  0.00025  seconds 
after  the  perturbation.  At  this  time,  the  largest  difference  in  Clrm  between  the  full  and 
partially-filled  cylinders  is  computed  (see  Figure  26).  The  sidewall  component  of  the  shear 
stress,  and  thus  Clrm,  is  proportional  to  the  radial  derivative  of  the  azimuthal  velocity 
( dv/dr  —  v/r).  The  computed  difference  in  this  derivative  for  the  two  cylinder-fill  ratios, 
when  integrated  over  the  entire  sidewall,  yields  the  computed  difference  in  the  values  of 
Clrm- 

Figure  29  shows  the  relative  difference  in  the  computed  Clrm  for  the  filled  and 
partially-filled  cylinders.  For  this  plot  the  computation  was  extended  to  2000  time  steps 
or  0.005  secs.  Whereas,  the  difference  in  Clrm  between  /  =  1  and  /  =  0.95  approaches 


zero  in  the  limit,  the  relative  difference  in  Clrm  approaches  a  constant  that  is  greater  than 
0.3  for  this  case.  Since  the  relative  difference  is  increasing  with  time,  Clrm  for  /  =  1  is 
approaching  zero  faster  than  the  difference  in  Clrm  for  /  =  1  and  /  =  0.95. 

VIII.  CONCLUSIONS 

The  internal  liquid  flowfield  in  a  partially-filled  right  circular  cylinder  in  rapid  axial 
rotation  has  been  investigated.  Currently,  Reynolds  numbers  up  to  100,000  have  been  com¬ 
puted.  The  present  work  uses  the  axisymmetric  Navier-Stokes  equations  but  is  restricted 
to  spin-up  of  the  liquid  from  one  state  of  solid  body  rotation  to  another.  Limits  on  the  ax¬ 
ial  rotation  rate  of  the  container  are  chosen  such  that  the  shape  of  the  free  surface  remains 
vertical.  The  present  work  uses  the  stream  function-vorticity  form  of  the  unsteady,  viscous, 
incompressible  Navier-Stokes  equations.  As  a  result,  the  finite-difference  analog  of  these 
equations  are  directly  solved,  using  Gauss-Seidel  SOR,  without  modification  with  artificial 
compressibility  terms.  Upwind  differencing  has  been  selectively  employed  for  the  convec¬ 
tive  terms  of  the  governing  equations.  Logarithmic  grid  stretching  has  been  employed  near 
wall  boundaries  of  the  flow,  as  well  as  the  free  surface. 

Verification  of  the  spatial  and  temporal  accuracy  of  the  numerical  scheme  has  been 
demonstrated  for  a  rotating  cylinder  full  of  liquid.  The  required  digital  computer  resources, 
such  as  run  time  and  storage,  are  minimal  with  typical  cases  averaging  30  CPU  minutes 
on  a  VAX  8600  mini-computer  to  attain  steady-state  solutions  (in  addition  to  the  time 
accurate  solutions  saved  during  the  run). 

Free  surface  boundary  conditions  on  vorticity  and  circulation  were  derived  based  on 
the  assumption  of  zero  tangential  shear  stress.  A  Dirichlet  boundary  condition  was  chosen 
for  the  stream  function  that  was  consistent  with  the  assumption  of  a  vertical  free  surface. 
These  conditions  were  also  employed  by  Shadday.9  Both  computations  show  an  inadequate 
correlation  with  LDV  measurements  of  the  axial  velocity  component  at  the  free  surface. 
The  present  work  employs  a  computational  scheme  based  on  the  stream  function-vorticity 
form  of  the  governing  equations  that  is  very  different  from  the  scheme  used  by  Shadday. 
In  addition,  the  present  work  uses  computational  grids  that  are  highly  clustered  at  the 
boundaries  of  the  domain.  The  similar  results  generated  by  these  two  very  different  com¬ 
putational  approaches  points  to  the  free  surface  boundary  conditions  as  the  likely  source 
of  discrepancy.  However,  a  possible  velocity  bias  in  the  LDV  measurements  near  the  free 
surface  was  described  by  Shadday. 

The  relative  importance  of  the  free  surface  was  investigated  by  comparison  of  a 
Reynolds  number  1000  flow  in  a  full  and  partially-filled  cylinder  during  spin-up  from  one 
state  of  solid  body  rotation  to  another.  The  sidewall  component  of  the  liquid  roll  moment 
for  the  partially-filled  cylinder  differed  significantly  from  that  of  the  cylinder  full  of  liquid. 
The  endwall  component  was  virtually  unchanged. 
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Figure  3.  Finite-difference  grid  (31  x  31)  with  logarithmic  stretching  (d  =  b  =  0.05)  for 
full  cylinder  ( c/a  —  1.0) 
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Figure  4.  Finite-difference  grid  (31  x  31)  with  logarithmic  stretching  (d  =  b  =  0.05)  for 
70%  filled  cylinder  ( c/a  =  1.0) 
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Figure  5.  Azimuthal  velocity  profiles  at  cylinder  midplane  for 
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Figure  9.  Inertial  oscillations  during  spin-up  from  a  previous  state  of  rigid  body  rotation 
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Figure  10.  Inertial  oscillations  during  spin-up  from  a  previous  state  of  rigid  body  rotation 
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Figure  11.  Inertial  oscillations  during  spin-up  from  a  previous  state  of  rigid  body  rotation 
at  r  =  0.75,  2  =  0.3182 
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Figure  13.  LDY  measurement  stations  for  liquid  centrifuge  experiment  (Ref.  9);  1:  z 
0.497,  2:  z  =  0.9.  3 :  z  =  1.304,  4:  z  =  1.706 
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Figure  14.  Finite-difference  grid  (30  x  61)  for  Re  =  93518,  a  =  1.028,  36.7%  filled. 
(z^~)endwalh  =  0.0018.  ( ^r)sideuall  =  (^r)free  surface  =  0.0014 
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Figure  17.  Axial  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  z 


Radial  Distribution  of  Azimuthal  Velocity 

Z/a  =  0.9 

1  o  LDV  Data  (Ref.  9) 

Com pu tat i on  (Re f .  9} 

Present  Computation  ^*5 


0.795  0.836  0.877  0.918  0.959 


Figure  18.  Azimuthal  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled, 
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Figure  19.  Axial  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled, 
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Figure  21.  Axial  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  2 


lh  uv\ -  v  wnvwnwwvmwvRwjwivi  F.J1  . 


Radial  Distribution  of  Azimuthal  Velocity 

Z/a  =  1.706 


1 


0.95 


0.90 


0.85 


0.80 


0.795  0.836  0.877  0.918  0.959 


r 


Figure  22.  Azimuthal  velocity  profile  for  Re  =  93518,  a  =  1.028,  36.7%  filled,  2 
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Figure  23.  Finite-difference  grid  (21  x  41)  for  Re  =  1000,  a  =  3,  full  cylinder 


Figure  24.  Finite-difference  grid  (21  x  41)  for  Re  =  1000,  a  =  3,  95%  filled  cylinder 
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Figure  29.  Relative  difference  in  total  liquid  roll  moment  coefficient  for  / 
0.95  vs  time  for  Re  =  1000.  a  =  3 
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List  of  Symbols 


a  cylinder  radius 

b  grid  stretching  parameter  for  the  z-direction  (see  Eq.  35) 
c  cylinder  half-height 

d  grid  stretching  parameter  for  the  r-direction  (see  Eq.  34) 

Clrm  coefficient  of  liquid  roll  moment 

E  Ekman  number,  v/ftjc2,  the  ratio  of  viscous  to  Coriolis  force. 

/  liquid  fill  ratio  of  the  cylinder,  Z,/(2c) 

Fr  Froude  number,  Cl2a/g 

g  gravitational  acceleration,  9.81  ft/sec2 

P, p  dimensional,  nondimensional  local  pressure,  inertial  reference  frame,  p  =  P/pfl^a2 

R,r  dimensional,  nondimensional  radial  coordinate,  r  =  R/a 

ri  radial  location  of  the  free  surface  at  z  =  0  (cylinder  bottom  wall) 

r2  radial  location  of  the  free  surface  at  z  =  2a  (cylinder  top  wall) 

Re  Reynolds  number,  Qja2/ v 

ri  radial  location  of  the  vertical  free  surface  (see  Eq.  6) 

T,  r  dimensional,  nondimensional  time,  r  =  TQf 

U, u  dimensional,  nondimensional  radial  velocity  component,  u  =  U/Clja 

V,  v  dimensional,  nondimensional  azimuthal  velocity  component,  v  =  V/f If  a 
W}  w  dimensional,  nondimensional  axial  velocity  component,  w  =  W/fi/a 

Z ,  z  dimensional,  nondimensional  axial  coordinate,  z  =  Z/a 
Zt  liquid  level  at  the  cylinder  spin  axis  for  zero  axial  spin. 

Z0  liquid  level  at  the  cylinder  spin  axis  for  Qj  spinrate. 

Greek  Symbols 

q  cylinder  aspect  ratio,  c/a 

T, '/  dimensional,  nondimensional  circulation,  7  =  T/Qja 2 

A r  difference  in  free  surface  radial  location  between  top  and  bottom  cylinder  walls,  r2  —  rj 

Z.C,  dimensional,  nondimensional  vorticity,  £  =  Z/Sl] 

9  azimuthal  coordinate 

p  dynamic  viscosity  of  the  liquid  fill 

v  kinematic  viscosity  of  the  liquid  fill 
p  density  of  the  liquid  fill 

T0t  tangential  shear  stress  in  the  0-direction  on  a  surface  with  normal  in  the  r-direction 

tzt  tangential  shear  stress  in  the  2-direction  on  a  surface  with  normal  in  the  r-direction 

ty,  v  dimensional,  nondimensional  stream  function,  ip  —  'i/Qja3 

Q,  initial  axial  spinrate  of  the  cylinder  (rad/sec) 

SI j  final  or  steady  state  axial  spinrate  of  the  cylinder  (rad/sec) 

Qj  axial  spinrate  of  the  top  wall  of  the  cylinder  (rad/sec) 
us  Gauss-Seidel  relaxation  parameter 


Subscripts 

r  first  derivative  with  respect  to  the  radial  coordinate 
rr  second  derivative  with  respect  to  the  radial  coordinate 
2  first  derivative  with  respect  to  the  axial  coordinate 
zz  second  derivative  with  respect  to  the  axial  coordinate 


Superscripts 
k  current  iteration 
k  -f  1  successive  iteration 
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